#################################################################
## Military aid, regime vulnerability, and political violence  ##
##                                                             ##
## Andrew Boutton                                              ##
##                                                             ##
## Appendix Figures 1 & 4                                      ##
##                                                             ##
## Version: 7/24/2019                                          ##
#################################################################

data = read.dta("bjps.dta")

## Figure 1 -- Regime type frequency over time ##

regime.yr <- ddply(data, .(year), summarise,
                   democracy = sum(democ, na.rm = TRUE),
                   party = sum(party, na.rm = TRUE),
                   personalist = sum(personal, na.rm = TRUE),
                   military = sum(milregime, na.rm = TRUE),
                   monarchy = sum(monarchy, na.rm = TRUE),
                   other = sum(other, na.rm = TRUE) )
b <- barplot(t(as.matrix(regime.yr[,2:7])), plot = FALSE) 
a <- b[seq(5, length(b), 5)] 
par(cex.axis = .60, mai = c(0.5, 0.25, 0.25, 0.1))
barplot(t(as.matrix(regime.yr[,2:7])), border = "black", space = 0.2,
     axes = FALSE, ylim = c(0,175),
     col = c("blue", "midnightblue", "dodgerblue", "skyblue", "turquoise4", "lightgreen"),
     legend = c("Democracy", "Dominant party", "Personalist", "Military regime", "Monarchy", "Other"),
     args.legend = list(x = "top", y = 150, bty = "n", border = "white", cex = 0.95))  
axis(1, at = a, labels = c(1950,1955,1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010), tick = TRUE, pos = -1, cex.axis=1)
axis(2, las = 1, tick = TRUE, pos = 0,cex.axis=1)     

## Figure 4 - Interaction plot 

interaction_plot_continuous <- function(model, effect, moderator, interaction, varcov="default", minimum="min", maximum="max", incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=80, rugplot=T, histogram=T, title="Marginal effects plot", xlabel="Value of moderator", ylabel="Estimated marginal coefficient"){

  # Define a function to make colors transparent
  makeTransparent<-function(someColor, alpha=alph){
    newColor<-col2rgb(someColor)
    apply(newColor, 2, function(curcoldata){rgb(red=curcoldata[1], green=curcoldata[2],
      blue=curcoldata[3],alpha=alpha, maxColorValue=255)})
  }

  # Extract Variance Covariance matrix
  if (varcov == "default"){
    covMat = vcov(model)
  }else{
    covMat = varcov
  }

  # Extract the data frame of the model
  mod_frame = model.frame(model)

  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]

  # Set range of the moderator variable
  # Minimum
  if (minimum == "min"){
    min_val = min(mod_frame[[moderator]])
  }else{
    min_val = minimum
  }
  # Maximum
  if (maximum == "max"){
    max_val = max(mod_frame[[moderator]])
  }else{
    max_val = maximum
  }

  # Check if minimum smaller than maximum
  if (min_val > max_val){
    stop("Error: Minimum moderator value greater than maximum value.")
  }

  # Determine intervals between values of the moderator
  if (incr == "default"){
    increment = (max_val - min_val)/(num_points - 1)
  }else{
    increment = incr
  }

  # Create list of moderator values at which marginal effect is evaluated
  x_2 <- seq(from=min_val, to=max_val, by=increment)

  # Compute marginal effects
  delta_1 = beta_1 + beta_3*x_2

  # Compute variances
  var_1 = covMat[effect,effect] + (x_2^2)*covMat[interaction, interaction] + 2*x_2*covMat[effect, interaction]

  # Standard errors
  se_1 = sqrt(var_1)

  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf)/2))
  upper_bound = delta_1 + z_score*se_1
  lower_bound = delta_1 - z_score*se_1

  # Determine the bounds of the graphing area
  max_y = max(upper_bound)
  min_y = min(lower_bound)

  # Make the histogram color
  hist_col = makeTransparent("grey")

  # Initialize plotting window
  plot(x=c(), y=c(), ylim=c(min_y, max_y), xlim=c(min_val, max_val), xlab=xlabel, ylab=ylabel, main=title)

  # Plot estimated effects
  lines(y=delta_1, x=x_2)
  lines(y=upper_bound, x=x_2, lty=2)
  lines(y=lower_bound, x=x_2, lty=2)

  # Add a dashed horizontal line for zero
  abline(h=0, lty=3)

  # Add a vertical line at the mean
  if (mean){
    abline(v = mean(mod_frame[[moderator]]), lty=2, col="red")
  }

  # Add a vertical line at the median
  if (median){
    abline(v = median(mod_frame[[moderator]]), lty=3, col="blue")
  }

  # Add Rug plot
  if (rugplot){
    rug(mod_frame[[moderator]])
  }
  # Add Histogram (Histogram only plots when minimum and maximum are the min/max of the moderator)
  if (histogram & minimum=="min" & maximum=="max"){
    par(new=T)
    hist(mod_frame[[moderator]], axes=F, xlab="", ylab="",main="", border=hist_col, col=hist_col)
  }
}

## Estimate count model
## 'dp2' calculates the annual relative coup threat. More formally, dp2 = annual coup probability - country mean of coup probability. Adapted from Koga 2016.

m4 <- glm.nb(attks ~  attkst_1 + ln_securityMA + dp2 + dp2Xsecurity + personalt_1 + democt_1 + milregimet_1  + partyt_1 +
                  monarchyt_1 + ln_gdpt_1 + ln_popt_1 + civwar + rivalryt_1 + mediascore + post911 + coldwar + as.factor(ccode) + as.factor(year), data=data)

## Plot
interaction_plot_continuous(m4, effect="ln_securityMA", moderator="dp2", interaction="dp2Xsecurity", mean=T,
                             title="Regime-change coups",xlabel="Relative coup threat",
                             ylabel="Marginal effect of military aid on violence")